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Solving two-dimensional coupled Burgers 


equations via a stable hybridized 
discontinuous Galerkin method 


S. Baharlouei, R. Mokhtari*© and N. Chegini 


Abstract 


The purpose of this paper is to design a fully discrete hybridized discon- 
tinuous Galerkin (HDG) method for solving a system of two-dimensional 
(2D) coupled Burgers equations over a specified spatial domain. The semi- 
discrete HDG method is designed for a nonlinear variational formulation 
on the spatial domain. By exploiting broken Sobolev approximation spaces 
in the HDG scheme, numerical fluxes are defined properly. It is shown that 
the proposed method is stable under specific mild conditions on the stabi- 
lization parameters to solve a well-posed (in the sense of energy method) 
2D coupled Burgers equations, which is imposed by Dirichlet boundary 
conditions. The fully discrete HDG scheme is designed by exploiting 
the Crank—Nicolson method for time discretization. Also, the Newton— 
Raphson method that has the order of at least two is nominated for solving 
the obtained nonlinear system of coupled Burgers equations over the rect- 
angular domain. To reduce the complexity of the proposed method and the 
size of the linear system, we exploit the Schur complement idea. Numerical 
results declare that the best possible rates of convergence are achieved for 
approximate solutions of the 2D coupled Burgers equations and their first- 
order derivatives. Moreover, the proposed HDG method is examined for 
two other types of systems, that is, a system with high Reynolds numbers 
and a system with an unavailable exact solution. The acceptable results 
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of examples show the flexibility of the proposed method in solving various 
problems. 


AMS subject classifications (2020): 65M60, 65M12 


Keywords: Coupled Burgers equations; hybridized discontinuous Galerkin 
method; stability analysis. 


1 Introduction 


Throughout the history of science, finding the analytical and especially nu- 
merical solutions of nonlinear evolution equations such as Burgers and cou- 
pled Burgers equations [3, 28, 31, 35], KdV type equations [2, 4, 27], Navier— 
Stokes equations [34], and nonlinear Schrédinger equations [7] play crucial 
roles in various fields of science and engineering for the detection of physical 
phenomena. The system of two-dimensional (2D) coupled Burgers equations, 
as a simplified form of some complex and practical equations in engineering 
such as the incompressible Navier-Stokes equation, is widely used in fluid 
dynamics such as modeling of the shock waves moving in viscous liquid [17], 
shallow water waves [18, 26], turbulent medium [5], and diffusion processes 
[1]. According to the new works that are done in the literature [24, 38, 39], 
we realize that providing methods of finding the numerical solutions of Burg- 
ers and coupled Burgers equations still have their importance. Moreover, 
numerical scientists consider Burgers and coupled Burgers equations as test 
problems to introduce and experiment with new numerical methods. In other 
words, these equations are used to compare different numerical methods in 
various aspects to choose and extend the most appropriate one to a spe- 
cialized subject. This paper proposes a stable scheme for solving the 2D 
nonlinear coupled Burgers equations over rectangular domains numerically. 


The general form of the 2D system of coupled Burgers equations reads as 


1 
uz + uu, + Duy — Ro (ee + Uyy) = 0, 


vb, + uv, + vv, (isa F Hips) = 0, 


Re 
or equivalently 


1 
u,-+U-Vu-— —Au=0, 
. (1) 


+ . — ij => 
ob. +U-Vo Re u= 0, 


where Re > 0 is the Reynolds number, U = (u,v)™, and x = (a,y) EN = 
(a,b) x (c,d) C R?. In this paper, system (1) is equipped by the Dirichlet 
boundary conditions and suitable initial conditions. 
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Analytical solution of system (1) can be obtained, for instance, by the 
Hopf-—Cole transformation; see [20]. Providing an explicit analytical solu- 
tion for system (1) is not trivial. However, if, by any chance, an explicit 
form becomes available, then evaluating the analytical solution requires high 
computational costs with a considerable amount of time, which may be ac- 
companied by uncontrollable errors regarding the discretization of the ana- 
lytical solution. Based on these reasons, it is requested design be stable and 
effective numerical methods for computing numerical solutions. For solving 
system (1), many numerical methods have been proposed, such as the de- 
composition method [19], Chebyshev spectral collocation method [23], and 
some others; for instance, see [24, 38, 39]. 

Since the main approach of this paper is directly related to the discon- 
tinuous Galerkin (DG) method and is considered a continuation of the local 
discontinuous Galerkin (LDG) method, it is necessary to briefly review the 
history and background of DG and LDG methods. The first DG method 
was proposed by Reed and Hill in 1973 for a time-independent linear hyper- 
bolic equation [6], and then it was utilized and developed for time-dependent 
partial differential equations (PDEs); see [10, 16]. Provable cell-entropy in- 
equality for L? stability, h-p adaptivity, and flexibility to handle complicated 
geometry for arbitrary order of accuracy with local in-data communication, 
and other abilities lead to applying the DG method to various types of differ- 
ential equations. To dominate the limitations of the DG method for solving 
high-order partial differential equations, an LDG method was proposed. This 
method was used for the first time for solving a second-order time-dependent 
convection-diffusion equation [15]. The main idea of the LDG method is the 
transformation of a high-order equation into a first-order system of equations 
before solving the new system by the DG method. Due to eliminating all of 
the auxiliary variables locally, the LDG method inherits all flexibilities of 
the DG method. Recent applications of the LDG method for higher-order 
nonlinear PDEs can be found, for instance, in [8, 25, 31]. 


The usage of the hybridization technique in the context of the finite ele- 
ment method goes back many years ago, while its application in the context 
of DG methods has a recent history and goes back to 2004. In fact, the 
hybridized discontinuous Galerkin (HDG) method was proposed for the first 
time by combining the DG method and continuous Galerkin (CG) method to 
solve the steady-state problems [11], and then it was generalized by Cockburn 
et al. [12, 13, 14]. Recently, HDG methods have been widely used to solve 
evolution equations numerically, in particular for compressible flow problems 
(22, 30, 33, 36, 37], Stokes flow [9, 21], continuum mechanics problems [29], 
and linear elasticity problems [32]. The HDG methods inherit the optimal 
convergence rate from the DG methods for approximate solutions and their 
derivatives with respect to spatial variables. HDG methods have two kinds 
of unknowns; global unknowns that are used in the definition of numerical 
traces (or in numerical fluxes) and obtained from the global system, and local 
unknowns that can be eliminated locally and are obtained by weak formula- 
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tion. Local and global unknowns are approximated by piecewise polynomials 
of degree k, respectively, in R? and R¢~!, where d is the dimension of the 
spatial domain. Due to the consideration of global unknowns, one can infer 
that the degree of freedom in the HDG method is reduced compared to the 
traditional implicit DG methods. The key to the success of the HDG method 
is the way of defining numerical fluxes that are based on global unknowns 
and stabilization parameters. The numerical fluxes of the HDG method are 
not defined uniquely in most situations, but those have to be defined in such 
a way that the desired definitions of numerical fluxes ensure the stability of 
the scheme. Also, the definitions of the numerical fluxes cause significantly 
smaller bandwidth than the corresponding matrices of the traditional CG 
method, and therefore lower computational cost is accessible in any HDG 
method. In solving a problem with nonsmooth solutions, the HDG method 
as a kind of DG method is a suitable scheme. This advantage is based on 
the fact that the HDG method produces numerical approximations using dis- 
continuous trial functions over the entire given domain. In summary, it is 
worth pointing out that the HDG method has unique properties, which make 
this method superior, such as reducing the degree of freedom compared to 
the traditional implicit DG methods, making smaller bandwidth compared to 
the corresponding matrices of traditional CG and DG methods, and having 
less computational time; see [4]. In this paper, we intend to use a kind of 
HDG method for discretizing the 2D coupled Burgers equations (1) in the 
spatial domain. 

The rest of the paper is organized as follows. In Section 2, some pre- 
requisites such as notations, discretization of temporal and spatial domains, 
and approximation spaces, are expressed in dimension two. Section 3 is ded- 
icated to the employment of the HDG method to the 2D coupled Burgers 
equations. In fact, in this section, a semi-discrete scheme is presented for the 
2D coupled Burgers equations with suitable definitions of numerical fluxes 
and stabilization parameters. In addition, the stability of the proposed semi- 
discrete HDG scheme is investigated in this section. In other words, we 
prove that the method is stable in the L? norm under certain conditions on 
the stabilization parameters. Then, a full discretization approach is designed 
in Section 4 by exploiting the Crank—Nicolson method for time discretization 
and Newton—Raphson as a nonlinear solver. Numerical experiments in Sec- 
tion 5 show that the optimal order of accuracy is derived by the proposed 
method. Also, by performing some experiments, the numerical solutions of 
system (1) are investigated for large Reynolds numbers. Moreover, a 2D 
problem with different values of Reynolds numbers is investigated such that 
its exact solution is unavailable. The conclusion is given in Section 6. The 
paper is ended with an Appendix. 
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2 Prerequisites 


Order to set up a system of weak formulation of coupled Burgers equations, it 
requires defining necessary notations and relevant approximation spaces for a 
desired HDG method. With T as a final time and for all t € (0, T], we consider 
a given bounded spatial domain 2 C R? with suitable partitioning. Suppose 
that the domain Q = (a,b) x (c,d) is split into conforming and uniform 
finite element meshes with N triangles such that in this mesh generation, all 
triangles have no intersection except in common edges or vertices. In general, 
each of these triangles is denoted by K. By considering h as the longest edge 
among triangles, the finite collection of disjoint elements, and the set of the 
boundaries of elements, respectively, are denoted by 


Hn ={K}, 0.4, = {OK}, 


where 2 = Ucey, K, and OK denotes the boundary of element K. The 
collection F;, = F? UF? is the set of all faces such that F? and F? represent, 
respectively, the set of interior and boundary faces. More precisely, the set 
of faces contains all edges of triangles. Let us consider two elements K~ and 
K+ and their common face e = OK~ NOK* € FP. As illustrated in Figure 1, 
n~ and n®™ are, respectively, the corresponding outward unit normal vectors 
of face e with respect to K~ and K+. Let v~ and vt be the limits of the 
function v at face e associated with OK+ and OK~, respectively. Thus the 
mean and jump values of an arbitrary real valued function v on the given 
face e are, respectively, defined as 


{oH} = 507 +o*), [el] = ea + otnt. 


We note that the mean and jump values of function v at boundary face 


Figure 1: Common face e of two elements K+,K~ with outward unit vectors. 
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e =0QNK € FP are determined as {{v}} = v and [[v]] = un, respectively. 
So, the mean and jump of function v can be rewritten as 

— f (ut +v7)/2,€€ FP, _ futnt+u-n,e€ FR, 

eS {" e€ FP, Mal: un, e€ FP. 


To obtain weak formulations of the 2D coupled Burgers equations, one needs 
to define appropriate approximation spaces. Regarding the nature of any 
DG method, broken Sobolev spaces are relevant spaces for approximating 
the solutions of system (1) via the HDG method. The corresponding broken 
Sobolev space, associated with the partition .%, is defined as 


H'(.%) = {v:Q9R:0 |ce H'(K), for all K € AH}, 
and associated with the set F;, is defined as 
M! (Fp) = {u: Fr 7 R: pw |e€ H1(e), for all e € Fy}. 


Discontinuous finite element spaces for scalar and vector valued functions, as 
subspaces of broken Sobolev space H1(.%,) are, respectively, defined by 


See {w © H(%,) : w |e Pa(K), for all K€ Hii, 
—_— {w € (H1(%))? : w lee (Pr(K))2, for all K€ Hi, 


where Px(K) is the set of polynomials of degree at most k on the element 
K € %,. The approximation space of the broken Sobolev space over Fp, (or 
skeleton space) is defined as 


Minn ={u€ M' (Fh) : uw |e€ Pr(e), for all e € Fy}. 


Regarding the boundary conditions, it is needed to define the appropriate 
subspace of the skeleton space. Consider Dirichlet boundary conditions and 
the boundary data b, and by on OQ, which are associated with u and », re- 
spectively. Let [, and I’, be collections of boundary faces in which boundary 
data b, and by are specified over [’,, and [',, respectively. Based on the given 
boundary conditions, we define 


Mix (l.T) = {yu € My, + w(x) = T(x), x eT}, 


where I € {Ty,Ty}, and I is the L? projection with respect to the skeleton 
space of the boundary of the domain 2. The approximation spaces Sp,%, Sn,x, 
and M),, are equipped by the following inner products, respectively, 


(wi, wa).%, = > (wi,wa)e, (in badow, = Y> (ua, He)an; 
KEHn KEK, 
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where 
(wi, Wa)ic =|} w1(x) - W2(x) dx, (1, M2)aK =| fli + 2 ds, 
K OK 


in which w,, we are defined on .%, and p41, ple are defined on 0.%,. By 
considering vector functions w = (wi, w2)', z = (21, 22)7, w= (t1, f2)T, and 
7 = (m,12)7, the following inner products are needed 


(w, Z) x), = x (w,Z)K, (L, 1) aX, = S- (4, N)aKk; 


KExX, KEK, 


where w 1, w2, 21, and zz are defined on .%, and 1, 2, M1, and 72 are 
defined on 0.%,. Besides, we have 


(w,z)K = ((wi, 21)K, (we, 22)x)", (Hs N)aK = ((H1,71)aK; (H2,2)aK)". 


3 Construction of the semi-discrete HDG method 


As mentioned, we assume that system (1) is equipped by the Dirichlet bound- 
ary conditions over the rectangular domain 2. The initial step is to refor- 
mulate the 2D coupled Burgers equations (1) into a first-order system of 
equations. By defining the auxiliary variables P = (p1,p2)™ = (Vu)™ and 
Q = (q1, 42)" = (Vv)', the corresponding first-order system of (1) reads as 


1 
u+U-Vu-—V-P=0, 
P V 0 
— us ; 
1 (2) 
vb, +U- Vo - —V-Q=0, 
Re 
Q-Vv=0. 


By establishing the corresponding semi-discrete HDG method of the system 
(2), the stability of the semi-discrete method over the temporal interval [0, ¢] 
for t € (0, T], is explained in the next subsection. 


To have a corresponding conditionally well-posed problem of the system 
(2), it is worth pointing out that this system should be equipped with initial 
and boundary conditions. Weak formulation of the system (2) can be formed 
by multiplying each equation of (2) by an appropriate test function, integrat- 
ing over each element K € .%,, and using the Green’s first identity. Conse- 
quently, the aim is to find numerical approximations (u,v, P,Q) € Si “x Si. k 
such that for all test functions (wi, w2, Wi, W2) € Sik x Sik and K € &, it 
holds that 
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(4, Wie + (U- Vu, wide + (GP Vwi)x + (-——Pn, w1)9 
(Pw) e+ ((u, Vor) ac — (an, w:))o 


1 
(v;, Wa) + (U- V0, we) + (Re@: Vw) + (“750m Ww2)a 


e Under imposed boundary conditions, numerical traces t € Mp7 (bu, Pu) 
and 6 € Mn.x(bo, Ty) are properly defined for all z € F;, as 


fe — Jou, zETy, 7 _ J bp, zET», 
we)={e rer, Ly, y= fe foe iy 


where (€,¢) € Mp. (0,0u) x Mn.x(0,P,) is a global unknown pair. It 
can be observed that boundary data b, and by, are imposed in the 
definitions of the numerical traces @ and @, respectively, on Py and Ty. 
One can infer that & and @ are global unknowns corresponding to the 
faces without a defined boundary data. 


e In order to guarantee the stability of the semi-discrete method, numer- 


a oe 
ical fluxes —p-P and —,-Q are defined as 


ie 1 . ie 1 : 
Re = TR + Tu wa, Rot = Ree + oly — O)n, 


(5) 
where n is the outward unit normal vector with respect to the considered 
face. In (5), 7 and o are the stabilization parameters. The valid range 
of parameters 7 and o are determined in the stability theorem of the 
proceeding subsection. We note that the definitions of the numerical 
fluxes in (5) are not unique and depend on the form and physics of the 
problems. 


Remark 1. It is noteworthy that numerical fluxes and stabilization parame- 
ters play a key role in the stability of the semi-discrete method. We emphasize 
that functions —_P and -#9 on each element edge are approximated by 
their corresponding numerical fluxes so that the numerical fluxes are single- 
valued continuous functions across the element edges. In HDG methods, 
numerical fluxes depend on the numerical traces while global unknowns in 
the definitions of numerical traces depend on the faces. 


Due to the fact that @ and @ contain two global unknown variables over 
(0, 7] x Q, two extra global equations on each face should be added to the 
system (3). The required global equations can be gained by enforcing the 
conservation of the fluxes. Thus, the global unknowns are obtained with the 
following extra global equations: 
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—- —- 


le? -nj]=0, for ec F®, eae -nj]=0, for ee F?. (6) 


Then the local unknowns u, v, P, and Q, can be found by solving weak 
formulation (3) in each element K € .%p,. 


3.1 Stability analysis 


In this subsection, we verify the numerical stability of the weak formulation 
(3) over the time interval (0, t], for all ¢ € (0, 7]. To do this, let homogeneous 
Dirichlet boundary conditions imposed to the weak formulation (3). We start 
the analysis by multiplying the first equation of (1) by u to get 


1 

= —y? 4 —U-V(u?) — Rota =0. (7) 
By integrating (7) over the given domain 2 and using the Green’s first iden- 
tity, we get 

1 1 1 Ou 

-— / wd 5 [Luv *) d a [vey ax- = [ —ds=0. (8 
ah PTS he eae a ia aed ron (8) 
By applying homogeneous Dirichlet boundary conditions to (8) and regarding 


| Vu- Vu dx > 0, one can conclude that (8) leads to 
Q 


ld 


1 
aire 2 al : 2 <0. 
7 Re mts EY V(u“) dx <0 (9) 


Integrating (9) over the time interval [0,t], for 0 < t < T, the following 
inequality holds: 


T 
ICO + f BH, dx <0), (10) 
with 
®@(t, A) = i U- V(t?) dx, 
A 
where t is the function of x and t, and A is a subdomain of 2. Also, Output 


of ®(t, A) is a function of variable t. In the same approach, from the second 
equation of (1), we get 


sa 
Ilo, lla +f B(v,0) dt < |v(-,0)|a. (11) 
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Theorem 1. Let weak formulation (3) be equipped by the homogeneous 
Dirichlet boundary conditions over the domain 2. By assuming 7 > 0 and 


a > 0, it can be proved that the solution of weak formulation (3) satisfies the 
following inequalities for all t € (0,7): 


T 

lu(-, Igy + | B(u, H,) dt < |ju(-,0)|[3g,, 
T 

lo) Ie +f B(v, Hh) at < lw, OI, 


Proof. By setting w, = u, wy = uP. w2 = v, and w2 = na, in the weak 
formulation (3) and summing the first three equations and the last three 
equations of (3) together, we get 


(12) 


ld nee | 

3 aul + SE IPIR: + Bie + $0(u,K) =0, 
1d ioe el 
sal’ Ile + Re! Cll Lok 5 

where 


Exc = ge(P, Vue + alu, VP) +(—geP M, U)aKk — Re (am, Plax, 
Fax = pg (Q, Vv) + 2G VQ)« + (59-0, v) aK — 5 (60, Q)ax. 


Using the divergence theorem, the following relations are obtained: 


RP, Vue + elu, VP)c = Re Ic V ) dz = = Ja-(Pu) -nds 
= e(P n ites (13) 
nm(Q, Vv) + mle, VQ) = Re ELV - (Qu) dz = 3 + Jax (Qu) -n ds 
= R(Q- n, UV) aK; 


By applying (13) to Fy and Fox, using 


(an, Pyar = (P ‘Nn, a)ax, (én, Q)ak = (Q ‘nh, O)aK; 
and adding 
Lp i)ak = 0 Le dak =0 
“\ Re N,U)ak = VY, Re N,v)ak =U, 


respectively, into E.x and E> x, we obtain 


Exvx = ere n,u— a)ak + (gP ‘nu- aja 

= (pe BP) -n,u- dor, 
(—#Q-n,v— d)ax + (RQ-n,v— B)aK 
( 


= ((-%Q+ BQ) -n,v—-d)ax. 
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Oe 


Using the definitions of —_P and —-9 from (5), we obtain 
Ey =<7,(u-@)? >ax, Ea =< 0,(u—6)? >ar . 
By summing F1,« and E2,x over all elements, we get 


Kex% E.x = Vi Kex, OTs (u _ a)? >axc=<T, (u — ti)? > aHys 
xen, £2.« = Uren, < % (v — 6)? >ax=< o,(v— 6) >ax, . 


According to the assumptions T > 0 and a > 0, we can conclude oicg Kp Fux > 
0 and ixex%, Eo > 0. Finally, by summing (12) over all elements, using 
the obtained results, and ||P||%y, , ||Q||3, => 0, we conclude 


d d 
Sllulee, + 8(u, Hi) $0, Slluly, + (0, i) <0. 


By integrating above relations over [0,¢] for all t € (0, 7], the assertion of the 
theorem is concluded. 


Remark 2. According to (10)—(11), by assuming 


a T 
i B(u,Q) dt >0, | B(v,Q) dt > 0, 
0 0 


one can verify that the 2D coupled Burgers equations (1) is well-posed in the 
sense of the energy method. Therefore, in this case and based on Theorem 
1, the proposed HDG method is stable with 7 > 0 and o > 0. 


Briefly, Theorem 1 and Remark 2 show that the proposed semi-discrete 
HDG method is stable for solving well-posed 2D coupled Burgers equations 
provided some specific mild conditions on the stabilization parameters. More- 
over, this stability is unconditional because we have no condition on the step 
sizes. 


4 Numerical algorithm and implementation issues 


In order to design a fully discrete approximation method for solving the 
2D nonlinear coupled Burgers equations (1), it is needed to apply a time- 
discretization approach to the weak formulation (3). To do this, we simply use 
the Crank—Nicolson method which is a method of order two. By considering 
time step At = 5 with J € N and time level t,, = nAt, for n = 0,..., J, the 
weak formulation (3) changes to 
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1 1 11 3. 

Aye wi + ae Vu", wide + a (ner Me + 9 (— RP) m wi)oKx = (wi), 
((P™,wi))c + (Cu, Vwi))c — ((a@"n,wi))ox = 0, 


i i A. al: a al 1 
(U"-Vo »wa)x + 5(p6@ »Vwa)e + 5((-R.@) n,we)ax = le(we), 


Re 
((Q”,w2))ic + ((v", Vwa))c — ((6"n, W2))ox = 0, 


1a 42 
Age eS 


where U" = (u",u")T, P® = (pt, pz)", and Q” = (qf, qf)", and 
ly(wi) = ay (u" 1, wi) + 5(U?* - Vu", wi) + (ge P", Vide 
13 ((—ReP)" "0, wi)ax, 


Io(we) = ar (v", wa) + 5(U"!- Vo", we) + §(q5Q” |, Vor) 


+3((— RQ)" *n, wa) ax. 


The superscripts n and n — 1 stand for the values at the time levels t¢, and 
tn—1, respectively. Likewise, the global equations should be considered at the 
time level ¢,,. By summing over all elements, inserting the flux definitions (5) 
into (6) and (14) at the time level t,, and also using boundary conditions (4), 
the algebraic system of equations or vector-matrix system can be obtained. 
The obtained system, steamed by exploiting the Crank—Nicolson method, 
is nonlinear, and we intend to solve it numerically so that preserves the 
second-order convergence in the temporal domain. Nevertheless, we exploit 
the Newton—Raphson method for solving the obtained nonlinear system. We 
set 


w” = (u™,u", pT Ph, at. gg €", 0") € SP. x Mp7 (0,0) x M),x(0,T»), 


where (w”,v”, py, po, W, G3, €", 6”) is the exact solution vector of system (14) 
and (6) at the time level t,,. With a suitable initial guess W,,,9, we are aiming 
to generate the following sequence of solution vectors 


Wri — Wri-l OW i i= 1, 2, cr) 


where W,,,; converges to the exact solution, namely, W”, as 2 tends to infinity. 
We note that 


OW = (OUn i; OUni, OPi nis OP2.n,i3 O91,n,i3 692, n,i3 dEn is OGn,i)s 


is obtained via the Newton—Raphson method. In the other words, 6W,,; 
is computed by solving the following linear variational formulation so that 
holds for all (w1,w2,w1,w2) € She x Sik and K € %, and (p1, p12) € 
Mh,x(0,T.) x Mh,x(0,T»): 
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41 (6Un;, W1) + a2(dUn 4, W1) + a3(dP1n,4,W1) + a4(dp2.ni,W1) 
+45 (06n,2, Wi) = (wi) 
1(0Un i; wii) ale b2(dP1,n,i> wii) Tr b3(d&n,i, wi) = I(wi 
4(OUn, és wi12) + bg (dp2,n,i5 w12) + bs (are wi2) = 15 (wie ’ 
C1(0Un i; we) + Co (5Un ji; we) + C3 (Od1,n,i; we) + C4 (52,n,is we) 
7 . +€5 (dEn, j,W2) = = a 
bi (dUn,i, wai) 4 + b2(dq1,n,i5 w21) ae b3(5Gn, as wa1) = ( 
ba(dun, io W22) 4 + bo (4g2,n,i5 W22) ah: bs (5Gn, io W22) = 6(W1 
+ d3 = ) = Ir( 
+ d3 ) ) = Is( 


lommon) 


1(0Un i, Ma) + do(6p1, nis 1) 4 (dp2,n,i5 LA rda(En, io M1 
1(0Un i; 2) + do(dqr, nis H2) 4 (d¢2,n,i; 2) — oda(SCn, iy M2 


d 
d 


where wy = (wi1,Wi2)', and we = (wa1,W22)'. To observe the definition 
of multilinear forms and linear functionals in (15), we refer the reader to 
Appendix of the paper. 


In order to solve the large and sparse linear variational formulation (15) 
more effectively, this system can be decomposed into two linear systems with 
smaller sizes by using the Schur complement idea. One can observe that (15) 
can be reformulated to the following vector-matrix equations: 

{ Miu Xn + Mi2Ynj = Ri, (16) 
M1 Xn + M22Yn i = Ra, 


where Xn i — [SUn i O0n i OP1 ni Opa,n,i bG1,n,i 6G2,n,i]", Yosi = [5En,i OGni]T 
are coefficients of approximate solutions with 


A; A. A3 Ay 0 0 As 0 
B, 0 B20 0 0 Bs 0 
B, 0 0B, 0 0 B, 0 
Mi=) a My = : 

11 Ci Cy, 0 0 Cs Ci , 12 0 Cs ’ 
0 B, 0 0 B 0 0 Bs 
0B, 0008B 0 By 

—rD4 0 TD, 0 Dp; Ds 0 
M. — oe M: —= ~ ~ ~ 
- | 0 os, | - "6 oD, 0 0 Dz Ds 

Ry = [Li Le Ls La Ls Le], Ro =[L7 Ls], 


In the above matrices and vectors, capital letters are interpreted as the ma- 
trix and vector representation of multi-linear forms and linear functionals 
defined in (15). Based on our experiences in the computer implementation of 
the HDG method, we have not benefited by not encountering non-invertible 
matrix M),. Regarding this fact, we assume that Mj), is invertible. Oth- 
erwise, it is not possible to propose the reduction of the complexity of the 
computations for solving (15), and so this system has to be solved directly. 
Based on the structure of the matrices in vector-matrix equations (16) and 
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the Schur complement issue, instead of solving (16), the following system of 
equations are solved in each iteration of the Newton—Raphson method: 


(Ma2 — Mo1Myy' Mi2)Yni = Re — Mar M7,'R1. (17) 
Thus X,,,, can be computed by 
Xn = MQ Ri — My MY. (18) 


Based on the Newton—Raphson approach, and Schur complement decom- 
position, we finish this section by representing the details of the designed 
HDG scheme in the following algorithm. 

Algorithm HDG algorithm for 2D coupled Burgers equations (1) 


Input: Spatial domain 2 and number of elements, namely, N, time inter- 
val [0,7] and number of time steps J, degree of approximate polynomials 
k, boundary data Ty and Ty, initial data, tolerance 0 < ¢, and stabilization 
parameters 7 and ao. 


Output: uw’, v7, pj, v3, a/, a, €7, and ¢7 that are the approximate solu- 
tions of ue. T); v(x, y,T), pif), po(a,y,T), qi(z,y,T), qa(z,y,T), 
(x,y, T) and ¢(2,y,T). 


Generate regular mesh for the domain Q. 
Set Wo by given initial and boundary conditions. 
Forn=1,2,...,J do 

Wno = Wn-1, )Wno = (e+1)1, i=0. 

While € < ||OW,, ;|| do 


Compute OW i+1 by Schur complement formulas (17) and (18) 


Writi = Writ OWnigi, i=it1 
end While 
W" = Wry 
end For 


5 Numerical results 


In this section, we aim to demonstrate the efficiency, validation, and ap- 
plicability of the proposed fully discrete HDG method for system (1). We 
observe that the semi-discrete HDG method for system (1) is stable over the 
time interval [0,¢], for all t € (0,T] provided that system (1) is well-posed in 
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the sense of energy method. To design a fully discrete version of the HDG 
method, we proposed an approach with the order of at least two for time 
discretization, that is Crank—Nicolson. Also, the Newton—Raphson method 
that has the order of at least two is proposed for solving the obtained non- 
linear system, and therefore, the loss of accuracy will not appear. As seen, 
to reduce the complexity of the proposed method and the size of the linear 
system, we exploited the Schur complement idea. Numerical experiments 
of the proposed HDG method are reported in three examples that they are 
selected from [35]. 

In Example 1, the 2D system (1) is considered to investigate the spatial 
order of accuracy of the proposed HDG method. Also, the results are reported 
for different Reynolds numbers. In Example 2, the HDG solution is examined 
for very high Reynolds numbers in the system (1). In Example 3, a 2D 
coupled Burgers equation without having any exact solution is solved. In 
this example, the HDG results are compared with the numerical results in 
[35] and [3]. 


Example 1. [35] Consider the 2D coupled Burgers equations (1) with Q = 
(0,1) x (0,1), T =1, and the following exact solutions: 


3 1 
"2 A+ exp(Be(—t — 4a + 4y)))’ 

3, 1 
4 4(1 + exp(22(-t — 4x + 4y))) 


The initial and boundary conditions can be derived from the exact solutions. 
In Table 1, L? error norms and corresponding orders are reported for Re = 1 
and tT = 0 = 0.5. As seen, satisfactory and high accuracy errors in Table 
1 indicate the good performance of our proposed method in solving system 
(1). Moreover, the results show the optimal convergence for approximate 
solutions u, v, and their first derivatives. As mentioned earlier, this optimal 
convergence is inherited from the DG method that is preserved well by our 
proposed method. In Table 2, the errors are reported for different Reynolds 
numbers Re = 0.1,1,10,100,200,500, approximate polynomials of degree 
k =2 and h = 0.2. For this test, we set 7 = o = 0.5 for Re = 0.1,1, 10, 100 
and +t = o = 2 for Re = 200,500. Note that, by increasing the Reynolds 
number, the effectiveness of dissipative terms in the system (1) will be elimi- 
nated gradually, and so we will face an inviscid system. Therefore, we expect 
that the accuracy of the method decreases as the Reynolds number increases. 
According to Table 2, we can observe that the proposed HDG method pro- 
duces acceptable approximate solutions even for high Reynolds numbers and 
the reduction of accuracy is acceptable. 

Here, we intend to do a test and check the dependence on the accuracy 
of the numerical solutions on the stability parameters. In Table 3, L? error 
norms and corresponding orders are reported for Re = 1 and tr = o = —0.5. 
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Table 1: L? error norms of approximate solutions wu, v, pi, p2, qi, and q2 together with 
their corresponding spatial orders of accuracy for Example 1 with Re = 1,7 =a =0.5 
at T= 1. 


h lu — ullo order | ||[pi —pilla order | |/p2—peallg order 
1 | 0.4 | 1.3059 E-6 3.7719 E-6 3.7620 E-6 

0.2 | 3.2810 E-7 1.99 | 9.3005E-7 2.02 | 9.3053 E-7 2.02 
0.1 | 7.6653 E-8 2.10 | 2.0728 E-7 2.017 | 2.0778 E-7 2.016 
2 | 0.4 | 1.2004 E-6 3.8744 E-6 3.8762 E-6 

0.2 | 1.5815 E-7 2.92 | 5.0946E-7 2.93 | 5.0959 E-7 2.93 
0.1 | 1.9734 E-8 3.00 | 6.3571E8 3.00 | 6.3588 E-8 3.00 
k h lv — dll order | ||qi1—qilla order | ||g2—qalla order 
1 | 0.4 | 1.3059 E-6 3.7719 E-6 3.7620 E-6 

0.2 | 3.2810 E-7 1.99 | 9.3005E-7 2.02 | 9.3053 E-7 2.02 
0.1 | 7.6653 E-8 2.10 | 2.0728 E-7 2.017 | 2.0778 E-7 2.016 
2 | 0.4 | 1.2004 E-6 3.8744 E-6 3.8762 E-6 

0.2 | 1.5815 E-7 2.92 | 5.0946E-7 2.93 | 5.0959 E-7 2.93 
0.1 | 1.9734 E-8 3.00 | 6.3571E-8 3.00 | 6.3588 E-8 3.00 


Table 2: L? error norms for Example 1 with approximate polynomial of degree k = 2 
and h = 0.1 for Re = 0.1,1, 10, 100, 250, and 500, at the final time T = 1. 


Re lu — ullor lpi —Pillo | llp2 — palle lv = dllo lar — gallo llaz — gallo 
0.1 1.6141 E-11 1.9859 E-10 1.9929 E-10 1.8628 E-11 2.7554 E-10 | 2.7644 E-10 
1 1.2378 E-8 3.9750 E-8 3.9774 E-8 1.2378 E-8 3.9751 E-8 3.9772 E-8 
10 6.7030 E-6 2.5580 E-5 2.6527 E-5 6.7030 E-6 2.5580 E-5 2.6527 E-5 
100 1.0638 E-3 1.9183 E-2 2.1982 E-2 1.0638 E-3 1.9183 E-2 2.1982 E-2 
200 3.3772 E-3 1.2515 E-1 1.3351 E-1 3.3773 E-3 1.2514 E-1 1.3351 E-1 
500 2.1209 E-2 7.1691 E-1 6.9343 E-1 2.1210 E-2 7.1691 E-1 6.9343 E-1 


We can observe that the HDG method with negative stabilization parameters 
produces numerical results with high and unacceptable errors. 


Example 2. [35] The aim of this example is to investigate the performance 
of the HDG method in solving system (1) with high Reynolds numbers. Con- 
sider system (1) with Q = (0,1) x (0,1), T = 1, and exact solutions 


=oat) cos(272) sin(y) 


Re(2 + exp(=32"*) sin(2mx) sin(y)) | 


2m exp( 


—5x7t 


2m exp(=Re 


Re(2 + exp(=32*) sin(27x) sin(y)) 


) sin(272) cos(7y) 


The errors of numerical solutions u and v are shown in Figures 2 and 3, 
respectively, for Re = 10000 and 100000. Note that, the results have been 
obtained by setting 7 = o = 20, h = 0.1, and k = 1. As mentioned in 
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Table 3: L? error norms of approximate solutions wu, v, pi, p2, q1, and q2 together with 


their corresponding spatial orders of accuracy for Example 1 with Re = 1, 7 =o = —0.5 
at T= 1. 
h lu — ullo order | ||[p1 —pilla order | |[p2—palla order 
1} 0.4 7.2644 2.4545 E+1 2.3701 E+1 
0.2 6.2532 0.22 | 1.4973 E+1 0.71 | 1.4567 E+1 0.70 
0.1 6.1053 0.03 | 1.4241 E+1 0.07 | 1.48304E+1 0.03 
2 | 0.4 | 2.9309 E+6 1.3503 E+5 1.6600 E+5 


0.2 | 8.4685 E+6 -1.53 | 8.4199 E4+5 = -2.6 8.0063 E+5 = -2.27 
0.1 | 2.2730 E+5 = -1.42 | 5.4620 E+4 = -2.70 | 4.6931 E+4 = -2.55 


k h lv — vllo order \la1 —qilla order llq2 —qa|lq order 
1 | 0.4 7.3625 2.6346 E+1 2.5501 E+1 
0.2 6.2372 0.24 1.4836 E+1 0.83 1.4870 E+1 0.78 
0.1 6.0686 0.04 1.4378 E+1 0.05 1.4228 E+1 0.06 
2 | 0.4 | 1.1331 E+6 5.3889 E+6 6.3966 E+6 


0.2 | 5.2046 E+6 = -2.20 | 5.1255 E+5 -3.25 | 4.9186 E+5 = -2.94 
0.1 | 9.5820 E+6 -0.88 | 2.3339 E+4 -2.19 | 1.9896 E+4 — -2.20 


Example 1, these high Reynolds numbers are going to omit the dissipative 
terms in system (1), but we can infer from Figures 2 and 3 that the behav- 
iors of approximate solutions still follow the exact solutions very well. This 
shows the flexibility and superiority of the proposed HDG method for solving 
different types of system (1) numerically. 


Example 3. [3, 35] In this example, a 2D problem with different values 
of Reynolds numbers will be investigated such that its exact solution is un- 
available. Consider the 2D coupled Burgers equations (1) over the domain 
Q = (0,0.5) x (0,0.5) with the initial conditions 


u(x, y, 0) = sin(aax) + cos(7x) o(z,y,0)=ax+y, 
and the boundary conditions 


(0,y,t) =cos(ry), u(0.5,y,t) = 1+ cos(ry), 
u(x,0,t) =1+sin(rz), u(#,0.5,t) = sin(7z) 
( t 


v(0,y,t)=y, v(0.5,y,t)=05+y, v(2,0,t)=2, v(2,0.5,t) =0.542. 


In the proposed HDG scheme, we set o = T = 2, h = 0.05, At = 0.001, and 
T = 0.625. According to this system that has no available exact solution, 
the only way to understand the correctness of the results is the compari- 
son them with the results of other papers. In the following, the results are 
compared with the results of [3, 35]. The numerical approximations u and 
v are illustrated in Figure 4 with k = 2 and Re = 50. Also, in Tables 4 
and 5, the numerical results are reported at some selected mesh points for 
Re = 50,500. We find that, the results of the proposed HDG method are in 
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Table 4: Comparison of computed values of u and v for Re = 50 for Example 3. Results 
are reported for approximate polynomials of degree two, h = 0.05, and At = 0.001 at 
final time T = 0.625. 


u v 
(x,y) HDG method [35] [3] HDG method [35] [3] 
0.1, 0.1) 0.96969 0.97146 0.96688 0.09817 0.09869 0.09824 
0.3, 0.1) 1.15072 1.15280 = 1.14827 0.14167 0.14158 0.14112 
0.2, 0.2) 0.86362 0.86307 0.85911 0.16915 0.16754 0.16681 
0.4, 0.2) 0.99136 0.97981 0.97637 0.18855 0.17109 0.17065 
0.1, 0.3) 0.66440 0.66316 0.66019 0.26491 0.26378 0.26261 
0.3, 0.3) 0.77587 0.77230 0.76932 0.24818 0.22654 0.22576 
0.2, 0.4) 0.59083 0.58180 0.57966 0.33124 0.32851 0.32745 
0.4, 0.4) 0.75273 0.75855 0.75678 0.38614 0.32499 0.32441 


Table 5: Comparison of computed values of u and v for Re = 500 for Example 3. Results 
are reported for approximate polynomials of degree two, h = 0.05, and At = 0.001 at 
final time T = 0.625. 


u v 
(x,y) HDG method [35] [3] HDG method [35] [3] 
(0.15, 0.1) 0.96114 0.96151 0.96650 0.08662 0.09230 0.09020 
(0.3, 0.1) 0.97324 1.03200 1.02970 0.07841 0.10728 0.10690 
(0.1, 0.2) 0.84445 0.87814 0.84449 0.17889 0.16816 0.17972 
(0.2, 0.2) 0.86926 1.06370 0.87631 0.16264 0.23690 0.16777 
(0.1, 0.3) 0.67883 0.67920 0.67809 0.26177 0.26268 0.26222 
(0.3, 0.3) 0.77557 0.79947 0.79792 0.21739 0.23550 0.23497 
(0.15, 0.4) 0.54874 0.58959 0.54601 0.31817 0.30419 0.31753 
(0.2, 0.4) 0.58850 0.78233 0.58874 0.30049 0.35294 0.30371 


good agreements with the presented results in [3, 35]. Hence, the proposed 
HDG method copes well with equations without the exact solution. 


6 Discussion and conclusion 


Numerical simulation of the 2D coupled Burgers equations via the HDG 
method has been studied in this paper, so that this system is equipped with 
appropriate initial and boundary conditions. In general, HDG methods have 
less computational time compared to the other DG methods, especially the 
LDG methods, which are the nearest to the HDG. The main reason for this 
advantage is the way of defining numerical fluxes. In the HDG method, the 
definition of numerical fluxes is not unique and depends on the form and 
physics of the problem. On the other hand, the stability of the method is 


Iran. j. numer. anal. optim., Vol. 13, No. 3, 2023,pp 397-425 


415 Solving two-dimensional coupled Burgers equations ... 


«10° 
3 
25 
15 
1 
os 


s 


; Error of approximate solution u(x, y, 1) 


baa? 


, Error of Sianeli solution u(x, y, 2) 


baee 


, Error of etctadeel solution v(x, y, 1) 


08 
0.6 

~ 
o4 
0.2 
ae 


Error of eseciadiieeata solution v(x, y, 2) 


1 

08 

0.6 
> 

0.4 

0.2 

aa 


ny 


a, 
ie 
Figure 2: The errors of approximate solutions u and v for Example 2 with Re = 10000 


at T = 1,2. The results are reported for 7 = o = 20, approximate polynomial of degree 
one, and h = 0.1. 
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Figure 3: The errors of approximate solutions u and v for Example 2 with Re = 100000 


at T = 1,2. The results are reported for tT = o = 20, approximate polynomial of degree 
one, and h = 0.1. 
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u(x, y, 0.625) 


v(x, y, 0.625) 


Figure 4: Approximate solutions u and v for Example 3 with Re = 50 at T = 0.625. 


The results are reported for 7 = o = 2, approximate polynomial of degree two with 
h = 0.05 and At = 0.001. 
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completely dependent on the fluxes. So, one of the hardships of using the 
HDG method is finding appropriate definitions of numerical fluxes that guar- 
antee stability. Fortunately, we presented a stable HDG method for solving 
system (1), while there is no stable (with a proven theorem) LDG method 
yet. Investigating the convergence of HDG methods for coupled and nonlin- 
ear problems is not easy. According to the authors’ knowledge, already, the 
convergence of the HDG methods has been studied just for some simple and 
linear equations. Convergence of the proposed method can be considered as 
one of our future works. 

The same as other HDG methods, by converting the initial system to a sys- 
tem of first-order equations and defining approximate broken Sobolev spaces 
associated with spatial partitioning, we set up the semi-discrete variational 
formulation of the coupled Burgers. Based on the structure of the HDG 
method, we have proposed numerical traces and fluxes for the variational 
formulation of the first-order equations. Numerical traces are supposed as 
global unknowns and depend on Dirichlet boundary data. Defining numer- 
ical flux in any HDG method plays a significant role in the stability of the 
semi-discrete HDG method over a time interval. After introducing appropri- 
ate numerical fluxes and imposing sufficient global equations over the spatial 
partitioning for system (1), the L? stability of the proposed semi-discrete 
HDG methods has been investigated under specific mild conditions on the 
stabilization parameters that are used in the definitions of numerical fluxes. 
With the intention of gaining a fully discrete scheme, the Crank—Nicolson 
method has been applied for time discretization. The choice of the Crank— 
Nicolson method was because of its unconditional stability and second-order 
accuracy. To preserve the second order of accuracy in time, the Newton— 
Raphson method has been nominated for solving the nonlinear system of 
equations. To solve the large and sparse linear variational systems, which 
is related to the Newton—Raphson method, the Schur complement idea has 
been used for reducing computational complexity and designing smaller sys- 
tems of equations. To explain the details of the HDG method, an algorithm 
has been prepared. To verify the efficiency of the proposed HDG method, the 
method was applied to some model problems. In the presented examples, we 
showed that approximate solutions and their first derivatives of degree k have 
converged at order k + 1, which is an optimal order of convergence. Also, in 
another example, the ability of the proposed HDG method was checked for 
solving the 2D coupled Burgers equations with different and high Reynolds 
numbers. Finally, we tested this method to solve a system without the ex- 
act solution, and pleasant results were observed. Regarding the flexibility of 
the method and numeric experiences, one can infer that the HDG method is 
one of the outstanding methods that has been exploited for various types of 
evolution problems in higher dimensions. 


Iran. j. numer. anal. optim., Vol. 13, No. 3, 2023,pp 397-425 


419 Solving two-dimensional coupled Burgers equations ... 


Appendix 


As mentioned in section 4, a fully discrete approximation method is obtained 
for solving the nonlinear coupled Burgers equations (1) by using the HDG 
and Crank—Nicolson methods for spatial and temporal discretization, respec- 
tively. Regarding the definitions of approximation spaces, the goal is to find 
(n,v”, pt, pS, ar,qe) € SS, and (€",¢") € Mn. (0,Tu) x Max (0,0), such 
that all equations in (6) and (14) are satisfied for n = 1,2,...,J. As stated 
in section 4, by summing over all elements, inserting the flux definitions (5) 
into (6) and (14) at the time level t,, and also using boundary conditions (4), 
the following system of equations is obtained: 


ice Wi) XH; ale 5 Up W1) th + 5 Uy, W1) x, — 2Re Far Pie te 1) Hh 
1 Ty. ie Tes 
~ Re app (Pay: W Wi) x, + 3 fu ,W1) ax, _ a6 W1) d.%t%\P un = 11 (w1), 
(pt, wi1)%, + (u”, (Wir)e)%,— < €'Ne, Wi >axy,\r, = l2(wi1), 
(py, wi2).4;, + (u”, (Wi2)y).%,— < E "Ny, Wi2 >a%,\r, = !3(wiz), 
i eee 1 Md, ys 
res , W2).%%, 15 2 (u” Up) W2)K a sv Uy) We), Re eran Cee 


o n nm 
~ Re apr (Gay: W2)x, + oe ,W2)a.%, — aC ,W2)a.x,\Py = la(we), 


(qf, wa1).0;, + (U", (war) ).%,— < OMe, Wai >a%,\Py = Is(war), 


(a3, Wa2).4;, + (UV, (we2)y).%,— < C "My, W22 >ax%,\r, = le(wa2), 


1 
Tu", Hr), \ FE ~ TE" Madar Fe — Re (Pi Be M1) ax, 


=F (PH My, M1)ax, = 0, 


1 
av", M2) at\F? — o(C”, M2) axty\F? a Ro Mt “Dy, 12) 0.x), 


—Re(@2 By, H2)ax%, = 9, 
(19) 
where n= (nz, ny)", 1 € Mp ,,(0), we © Mp,,(0), and 
1 n- I n-1, n— 1 n-1,n— 
li(w) = Ay * Ww) x, 3 (u : x£ * WwW) Hh, a a’ i * W) Xs 
iL ie 1 = T on, 
+ oR, (Piz *, wi), + We (Pu *,w1)#%, — ght 1 W)axh 
aoe 
+56 *, wax, + =(b, waxnr, 
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La(w) = 5500" wae, — (eT wan, — (0 HOP, w) an 
tclais wan + clay tw) — Flo", whoon 
+5 (O"* wi)oxm, + 5 (Bos W)axint, 

In(w) = (Dime, W)ax,rPur — Is(w) = (bimy, W)ax.nPws 


Is(w) = (bp Ne, W)a.%,rrysl6(w) = (dy Ny, W)ax,nPy- 


420 


Regarding to the nonlinear weak formulation (19), it is needed to convert 
this weak form to a linear variational form by a suitable iterative method. 
As explained in section 4, by applying the Newton—Raphson method to the 


nonlinear variational formulation (19), we intend to find 


OWn i _ (Uni; OUn ji; OP1njis Op2,n,is O41 ,n,i3 692,n,i3 dEn,is 5Gn,i)s 


such that for all (wi, W2, W141, W12, Wai; W22) E SP k and (41, [12) E Mp, (0, Pie) x 
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